#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>

#define NX 10
#define NY 10
#define NZ 10
#define nX NX+1
#define nY NY+1
#define nZ NZ+1
#define D 100
#define XM  NX*D
#define YM  NY*D
#define ZM  NZ*D

void ini(FILE *fp);
void dr_line(FILE *fp, float p[][3], int n, float *lc);
void dr_grid(FILE *fp, float *lc);
void dr_scale(FILE *fp, float *lc);
void dr_text(FILE *fp, float *tc, float *p);
void Viewpoint(FILE *fp);
void cal_pc_sub(float th, float pt[][3], float *vt, float *pcs, int *pcts);
void cal_pc(float th, float p[][3], float *v, float pc[][3], int *pct);
void dr_face(FILE *fp, float p[][3], int n, float *fc, float tr);
void mk_index_tri(int indt[][3]);
void mk_index_squ(int inds[][5][4]);
void mk_triang(FILE *fp, float pc[][3], int *pct, float *fc, float tr);
void mk_square(FILE *fp, float pc[][3], int *pct, float *fc, float tr);
void mk_index_face(int indf[][4]);
void mk_faces(FILE *fp, float pc[][3], int *pct, float *fc, float tr);
void count_pc(int *pct, int *num);
void search_f(int indf[][4], int self, int next, int *face, int *ptr);
void search_p(float pc[][3], int *pct, int indf[][4], int face, int ptr, float *pft, int *next);
void mk_equ_face_sub(FILE *fp, float th, float p[][3], float *v, float *fc, float tr);
void mk_cube(float x, float y, float z, float p[][3], float v[][nY][nZ], float *vt);
void mk_equ_face(FILE *fp, float th, float v[][nY][nZ], float *fc, float tr);
void filename(char *file);
void light(float *fc, float *lig);
void mk_data(float v[][nY][nZ]);
void cal_pc20(float th, float p[][3], float *v, float pc[][3], int *pct);
void cal_pc21(float th, float p[][3], float *v, float pc[][3], int *pct);
void count_pc2(int *pct, int *num);
void mk_faces2(FILE *fp, float pcs[][3], int *pcts, float *fc, float tr);
void cal_pc3(float th, float p[][3], float *v, float pc[][3], int *pct);
void mk_pt_index(int pind[][6]);
void mk_equ_face2_sub(FILE *fp, float th, float p[][3], float *v, float *fc, float tr);
void mk_equ_face2(FILE *fp, float th, float v[][nY][nZ], float *fc, float tr);
void clock_start(clock_t *start);
void clock_end(clock_t start);
void trans_v(float *v, float *vtr, float p[][3], float ptr[][3]);
void mk_cube2(float x, float y, float z, float v[][nY][nZ], float vt[][3][3]);
void mk_cube3(float x, float y, float z, float dx, float dy, float dz, float p[][3], float v[][3][3], float *vt);

void mk_equ_face001(FILE *fp, float th, float v[][nY][nZ], float *fc, float tr);
void rm_cube(float th, float *vt, int *rm);
void mk_equ_face0021(FILE *fp, float x, float y, float z, float th, float v[][nY][nZ], float *fc, float tr);
void mk_equ_face002(FILE *fp, float th, float v[][nY][nZ], float *fc, float tr);


void ini(FILE *fp)
{
	fprintf(fp, "#VRML V2.0 utf8\n");
}


void dr_line(FILE *fp, float p[][3], int n, float *lc)
{
	int i;

	fprintf(fp, "\nShape{");
	fprintf(fp, "appearance Appearance{");
	fprintf(fp, "material Material{");
	fprintf(fp, "diffuseColor 0 0 0");
	fprintf(fp, "}");
	fprintf(fp, "}");
	fprintf(fp, "geometry IndexedLineSet{");
	fprintf(fp, "coord Coordinate{");
	fprintf(fp, "point[");
	for(i=0; i < n ;i++) fprintf(fp, "%.0f %.0f %.0f ",p[i][0], p[i][1], p[i][2]);
	fprintf(fp, "]");
	fprintf(fp, "}");
    fprintf(fp, "coordIndex[");
	for(i=0; i < n ;i++) fprintf(fp, "%d ", i);
	fprintf(fp, "]");
    fprintf(fp, "color Color{");
	fprintf(fp, "color[");
    for(i=0; i < n; i++) fprintf(fp, "%f %f %f ",lc[0],lc[1],lc[2]);
	fprintf(fp, "]");
	fprintf(fp, "}");
    fprintf(fp, "colorIndex[");
	for(i=0; i < n ;i++) fprintf(fp, "%d ", i);
	fprintf(fp, "]");
	fprintf(fp, "}");
	fprintf(fp, "}");
}

void dr_grid(FILE *fp, float *lc)
{
	int j,k;
	float gp[2][3];

	for(j=0; j<NY+1 ;j++){
		for(k=0; k<NZ+1 ;k++){
			gp[0][0] = 0;
			gp[0][1] = j*(float)D;
			gp[0][2] = k*(float)D;
			gp[1][0] = (float)NX*(float)D;
			gp[1][1] = j*(float)D;
			gp[1][2] = k*(float)D;
			dr_line(fp, gp, 2, lc);
		}
	}
	for(j=0; j<NX+1 ;j++){
		for(k=0; k<NZ+1 ;k++){
			gp[0][0] = j*(float)D;
			gp[0][1] = 0;
			gp[0][2] = k*(float)D;
			gp[1][0] = j*(float)D;
			gp[1][1] = (float)NY*(float)D;
			gp[1][2] = k*(float)D;
			dr_line(fp, gp, 2, lc);
		}
	}
	for(j=0; j<NX+1 ;j++){
		for(k=0; k<NY+1 ;k++){
			gp[0][0] = j*(float)D;
			gp[0][1] = k*(float)D;
			gp[0][2] = 0;
			gp[1][0] = j*(float)D;
			gp[1][1] = k*(float)D;
			gp[1][2] = (float)NZ*(float)D;
			dr_line(fp, gp, 2, lc);
		}
	}
}

void dr_scale(FILE *fp, float *lc)
{
	int i,j,k,l;
	float sp[2][3];
	float p[3];

	for(j=0; j<nY ;j+=NY){
		for(k=0; k<nZ ;k+=NZ){
			sp[0][0] = 0;
			sp[0][1] = j*(float)D;
			sp[0][2] = k*(float)D;
			sp[1][0] = (float)NX*(float)D;
			sp[1][1] = j*(float)D;
			sp[1][2] = k*(float)D;
			dr_line(fp, sp, 2, lc);
			for(i=0; i<2; i++){
				for(l=0; l<3; l++){
					p[l] = sp[i][l];
				    dr_text(fp, lc, p);
				}
			}
		}
	}
	for(j=0; j<nX ;j+=NX){
		for(k=0; k<nZ ;k+=NZ){
			sp[0][0] = j*(float)D;
			sp[0][1] = 0;
			sp[0][2] = k*(float)D;
			sp[1][0] = j*(float)D;
			sp[1][1] = (float)NY*(float)D;
			sp[1][2] = k*(float)D;
			dr_line(fp, sp, 2, lc);
		}
	}
	for(j=0; j<nX ;j+=NX){
		for(k=0; k<nY ;k+=NY){
			sp[0][0] = j*(float)D;
			sp[0][1] = k*(float)D;
			sp[0][2] = 0;
			sp[1][0] = j*(float)D;
			sp[1][1] = k*(float)D;
			sp[1][2] = (float)NZ*(float)D;
			dr_line(fp, sp, 2, lc);
		}
	}
}

void dr_text(FILE *fp, float *tc, float *p)
{
	char s ='"';

	fprintf(fp, "\nTransform{");
	fprintf(fp, "translation %.0f %.0f %.0f ", (float)(p[0]-1.5*(XM/10)), p[1], p[2]);
	fprintf(fp, "rotation 1 0 0 1.57 ");
	fprintf(fp, "children[");
	fprintf(fp, "Shape{");
	fprintf(fp, "appearance Appearance{");
	fprintf(fp, "material Material{");
	fprintf(fp, "diffuseColor %f %f %f", tc[0], tc[1], tc[2]);
	fprintf(fp, "}");
	fprintf(fp, "}");
	fprintf(fp, "geometry Text{");
	fprintf(fp, "string[%c{%.0f %.0f %.0f}%c]", s, (p[0]/D), (p[1]/D), (p[2]/D), s);
	fprintf(fp, "fontStyle FontStyle{size %.0f}", (float)(XM/10));
	fprintf(fp, "}");
	fprintf(fp, "}");
	fprintf(fp, "]");
	fprintf(fp, "}");

}

void Viewpoint(FILE *fp)
{
	char s ='"';
	fprintf(fp, "\nViewpoint{position %.0f %.0f %.0f ", (float)(0.5*XM), (float)(-1.5*YM), (float)(0.5*ZM));
	fprintf(fp, "orientation 1 0 0 1.57 ");
	fprintf(fp, "description %centry%c}", s, s);
}

void cal_pc_sub(float th, float pt[][3], float *vt, float *pcs, int *pcts)
{
	int i;
	float a,b,s;

	if( ( (vt[0]>=th)&&(vt[1]>=th) )||( (vt[0]<th)&&(vt[1]<th) ) ){
		(*pcts) = 0; 
		for(i=0; i<3; i++){
			pcs[i] = 0;
		}
	}else{
		(*pcts) = 1;   /*点pc[]は存在する*/
		if( vt[0]==th ){
			for(i=0; i<3; i++){
				pcs[i] = pt[0][i];
			}
			return;
		}else if( vt[1]==th ){
			for(i=0; i<3; i++){
				pcs[i] = pt[1][i];
			}
			return;
		}else{
			s = vt[1]-vt[0];
	        a = th-vt[0];
	        b = s-a;
		    for(i=0; i<3; i++){
				pcs[i] = (b/s)*pt[0][i]+(a/s)*pt[1][i];
			}
		
		}
	}

}

void cal_pc(float th, float p[][3], float *v, float pc[][3], int *pct)
{
	float pt[2][3];
	float vt[2],pcs[3];
	int pcts;
	int i,j,k;

	for(k=0; k<2; k++){
		for(j=0; j<3; j++){
			for(i=0; i<3; i++){
				pt[0][i] = p[j+(k*4)][i];
		        pt[1][i] = p[(j+1)+(k*4)][i];
			}
	        vt[0] = v[j+(k*4)];
	        vt[1] = v[(j+1)+(k*4)];
	        cal_pc_sub(th, pt, vt, pcs, &pcts);
	        pct[j+(k*4)] = pcts;
	        for(i=0; i<3; i++){
				pc[j+(k*4)][i] = pcs[i];
			}
		}
	}

	for(j=0; j<2; j++){
		for(i=0; i<3; i++){
			pt[0][i] = p[(j*4)][i];
		    pt[1][i] = p[(j*4)+3][i];
		}
	    vt[0] = v[(j*4)];
	    vt[1] = v[(j*4)+3];
	    cal_pc_sub(th, pt, vt, pcs, &pcts);
	    pct[(j*4)+3] = pcts;
	    for(i=0; i<3; i++){
			pc[(j*4)+3][i] = pcs[i];
		}
	}

	for(j=0; j<4; j++){
		for(i=0; i<3; i++){
			pt[0][i] = p[j][i];
		    pt[1][i] = p[j+4][i];
		}
	    vt[0] = v[j];
	    vt[1] = v[j+4];
	    cal_pc_sub(th, pt, vt, pcs, &pcts);
	    pct[j+8] = pcts;
	    for(i=0; i<3; i++){
			pc[j+8][i] = pcs[i];
		}
	}

}

void dr_face(FILE *fp, float p[][3], int n, float *fc, float tr)
{
	int i;

	fprintf(fp, "\nShape{");
	fprintf(fp, "appearance Appearance{");
	fprintf(fp, "material Material{");
	fprintf(fp, "diffuseColor ");
	fprintf(fp, "%f %f %f",fc[0], fc[1], fc[2]);
	fprintf(fp, " transparency %.2f", tr);
	fprintf(fp, "}");
	fprintf(fp, "}");
	fprintf(fp, "geometry IndexedFaceSet{");
	fprintf(fp, "solid FALSE");
	fprintf(fp, " coord Coordinate{");
	fprintf(fp, "point[");
	for(i=0; i < n; i++) fprintf(fp, "%.0f %.0f %.0f ",p[i][0], p[i][1], p[i][2]);
	fprintf(fp, "]");
	fprintf(fp, "}");
	fprintf(fp, "coordIndex[");
	for(i=0; i < n;i++) fprintf(fp, "%d ",i);
	fprintf(fp, "]");
	fprintf(fp, "}");
	fprintf(fp, "}");
}

void mk_index_tri(int indt[][3])
{
	indt[0][0] = 0;
	indt[0][1] = 3;
	indt[0][2] = 8;

	indt[1][0] = 1;
	indt[1][1] = 0;
	indt[1][2] = 9;

	indt[2][0] = 2;
	indt[2][1] = 1;
	indt[2][2] = 10;

	indt[3][0] = 3;
	indt[3][1] = 2;
	indt[3][2] = 11;

	indt[4][0] = 4;
	indt[4][1] = 7;
	indt[4][2] = 8;

	indt[5][0] = 5;
	indt[5][1] = 4;
	indt[5][2] = 9;

	indt[6][0] = 6;
	indt[6][1] = 5;
	indt[6][2] = 10;

	indt[7][0] = 7;
	indt[7][1] = 6;
	indt[7][2] = 11;
}

void mk_index_squ(int inds[][5][4])
{
	inds[0][0][0] = 0;
	inds[0][0][1] = 2;
	inds[0][0][2] = 6;
	inds[0][0][3] = 4;

	inds[1][0][0] = 8;
	inds[1][0][1] = 9;
	inds[1][0][2] = 10;
	inds[1][0][3] = 11;

	inds[2][0][0] = 1;
	inds[2][0][1] = 3;
	inds[2][0][2] = 7;
	inds[2][0][3] = 5;

	inds[0][1][0] = 0;
	inds[0][1][1] = 1;
	inds[0][1][2] = 5;
	inds[0][1][3] = 4;

	inds[0][2][0] = 1;
	inds[0][2][1] = 2;
	inds[0][2][2] = 6;
	inds[0][2][3] = 5;

	inds[0][3][0] = 2;
	inds[0][3][1] = 3;
	inds[0][3][2] = 7;
	inds[0][3][3] = 6;

	inds[0][4][0] = 3;
	inds[0][4][1] = 0;
	inds[0][4][2] = 4;
	inds[0][4][3] = 7;

	inds[1][1][0] = 4;
	inds[1][1][1] = 6;
	inds[1][1][2] = 10;
	inds[1][1][3] = 9;

	inds[1][2][0] = 0;
	inds[1][2][1] = 2;
	inds[1][2][2] = 10;
	inds[1][2][3] = 9;

	inds[1][3][0] = 0;
	inds[1][3][1] = 2;
	inds[1][3][2] = 11;
	inds[1][3][3] = 8;

	inds[1][4][0] = 4;
	inds[1][4][1] = 6;
	inds[1][4][2] = 11;
	inds[1][4][3] = 8;

	inds[2][1][0] = 1;
	inds[2][1][1] = 3;
	inds[2][1][2] = 8;
	inds[2][1][3] = 9;

	inds[2][2][0] = 1;
	inds[2][2][1] = 3;
	inds[2][2][2] = 11;
	inds[2][2][3] = 10;

	inds[2][3][0] = 5;
	inds[2][3][1] = 7;
	inds[2][3][2] = 11;
	inds[2][3][3] = 10;

	inds[2][4][0] = 5;
	inds[2][4][1] = 7;
	inds[2][4][2] = 8;
	inds[2][4][3] = 9;

}

void mk_triang(FILE *fp, float pc[][3], int *pct, float *fc, float tr)
{
	int i,j,k;
	int indt[8][3];
	float ptri[3][3];
	int indtt[16];

	for(i=0; i<8; i++) indtt[i] = 0;
	mk_index_tri(indt);
	for(i=0; i<8; i++){
		if((indtt[i]==0) && ((pct[(indt[i][0])]==1) && (pct[(indt[i][1])]==1) && (pct[(indt[i][2])]==1))){
			for(k=0; k<3; k++){
				for(j=0; j<3; j++){
					ptri[k][j] = pc[(indt[i][k])][j];
				}
			    pct[(indt[i][k])] = 0;
			}
		    dr_face(fp, ptri, 3, fc, tr);
			indtt[i+4] = 1;
			i++;
		}
	}

}

void mk_square(FILE *fp, float pc[][3], int *pct, float *fc, float tr)
{
	int i,j,k,m,n;
	int inds[3][5][4];
	float psqu[4][3];

	mk_index_squ(inds);

	for(i=0; i<3; i++){
		for(j=1; j<5; j++){
		    if((pct[(inds[i][j][0])]==1) && (pct[(inds[i][j][1])]==1) && (pct[(inds[i][j][2])]==1) && (pct[(inds[i][j][3])]==1)){
			    for(k=0; k<4; k++){
				    for(m=0; m<3; m++){
						psqu[k][m] = pc[(inds[i][j][k])][m];
					}
			        pct[(inds[i][j][k])] = 0;
				}
		        dr_face(fp, psqu, 4, fc, tr);
				j++;
			}
		}
	}

	for(i=0; i<3; i++){
		n=0;
		for(m=0; m<4; m++){
			if(pct[(inds[i][0][m])]==1) n++;
		}
		if(n==4){
		        for(k=0; k<4; k++){
				    for(j=0; j<3; j++){
						psqu[k][j] = pc[(inds[i][0][k])][j];
					}
			        pct[(inds[i][0][k])] = 0;
				}
		        dr_face(fp, psqu, 4, fc, tr);
			    return;
		}
	}
}

void count_pc(int *pct, int *num)
{
	int i;

	(*num) = 0;
	for(i=0; i<12; i++){
		(*num)+=pct[i];
	}
}

void mk_index_face(int indf[][4])
{
	indf[0][0] = 0;
	indf[0][1] = 1;
	indf[0][2] = 2;
	indf[0][3] = 3;

	indf[1][0] = 4;
	indf[1][1] = 5;
	indf[1][2] = 6;
	indf[1][3] = 7;

	indf[2][0] = 0;
	indf[2][1] = 4;
	indf[2][2] = 8;
	indf[2][3] = 9;

	indf[3][0] = 1;
	indf[3][1] = 5;
	indf[3][2] = 9;
	indf[3][3] = 10;

	indf[4][0] = 2;
	indf[4][1] = 6;
	indf[4][2] = 10;
	indf[4][3] = 11;

	indf[5][0] = 3;
	indf[5][1] = 7;
	indf[5][2] = 8;
	indf[5][3] = 11;

}

void search_f(int indf[][4], int self, int next, int *face, int *ptr)
{
	int i,j;

	for(i=0; i<6; i++){
		if(i==self){
			continue;
		}
		for(j=0; j<4; j++){
			if(next==indf[i][j]){
				(*face) = i;
				(*ptr) = j;
			}
		}
	}

}

void search_p(float pc[][3], int *pct, int indf[][4], int face, int ptr, float *pft, int *next)
{
	int i,j;

	for(i=0; i<4; i++){
		if(i==ptr) continue;
			if(pct[(indf[face][i])]==1){
			    for(j=0; j<3; j++){
					pft[j] = pc[(indf[face][i])][j];
				}
		        (*next) = indf[face][i];
			}
	}
}

void mk_faces(FILE *fp, float pc[][3], int *pct, float *fc, float tr)
{
	float pf[6][3];
	float pft[3];
	int indf[6][4];
	int n;
	int i,j,k;
	int next;
	int face;
	int ptr;
	int self;

	count_pc(pct, &n);
	if(n<5) return;
	mk_index_face(indf);

	k = 0;
	self = 0;
	for(i=0; i<4; i++){
		if(pct[(indf[self][i])]==1){
			for(j=0; j<3; j++){
				pf[k][j] = pc[(indf[self][i])][j];
			}
			next = indf[self][i];
			k++;
		}
	}
	if(k==0){
		self = 1;
		for(i=0; i<4; i++){
		    if(pct[(indf[self][i])]==1){
			    for(j=0; j<3; j++){
					pf[k][j] = pc[(indf[self][i])][j];
				}
				next = indf[self][i];
			    k++;
			}
		}
	}

	for(k=2; k<n; k++){
		search_f(indf, self, next, &face, &ptr);
		search_p(pc, pct, indf, face, ptr, pft, &next);
		self = face;
		for(j=0; j<3; j++){
			pf[k][j] = pft[j];
		}
	}
	dr_face(fp, pf, n, fc, tr);

}

void mk_equ_face_sub(FILE *fp, float th, float p[][3], float *v, float *fc, float tr)
{
	float pc[12][3];
	int pct[12];



	cal_pc(th, p, v, pc, pct);

    mk_triang(fp, pc, pct, fc, tr);
    mk_square(fp, pc, pct, fc, tr);
	mk_faces(fp, pc, pct, fc, tr);   
}

void mk_cube(float x, float y, float z, float p[][3], float v[][nY][nZ], float *vt)
{
	int i;
	int xt,yt,zt;

	for(i=0; i<2; i++){
		p[0+(i*4)][0] = x;
		p[0+(i*4)][1] = y;
		p[0+(i*4)][2] = z+i*D;

		p[1+(i*4)][0] = x+D;
		p[1+(i*4)][1] = y;
		p[1+(i*4)][2] = z+i*D;

		p[2+(i*4)][0] = x+D;
		p[2+(i*4)][1] = y+D;
		p[2+(i*4)][2] = z+i*D;

		p[3+(i*4)][0] = x;
		p[3+(i*4)][1] = y+D;
		p[3+(i*4)][2] = z+i*D;
	}
	for(i=0; i<8; i++){
		xt = (int)p[i][0]/(int)D;
		yt = (int)p[i][1]/(int)D;
		zt = (int)p[i][2]/(int)D;
		vt[i] = v[xt][yt][zt];
	}

}

void mk_equ_face(FILE *fp, float th, float v[][nY][nZ], float *fc, float tr)
{
	float x,y,z;
	float vt[8];
	float p[8][3];

	for(x=0; x<(NX)*D; x+=D){
		for(y=0; y<(NY)*D; y+=D){
			for(z=0; z<(NZ)*D; z+=D){
				mk_cube(x, y, z, p, v, vt);
				mk_equ_face_sub(fp, th, p, vt, fc, tr);
			}
		}
	}
}

void filename(char *file)
{
	printf("output file: ");
	scanf("%s", file);
	getchar();

}

void mk_data(float v[][nY][nZ]){

	int i,j,k;
	float x,y,z;

	for(i=0; i<nX; i++){
		for(j=0; j<nY; j++){
			for(k=0; k<nZ; k++){
				x = i*(float)D;
				y = j*(float)D;
				z = k*(float)D;
				v[i][j][k] = (float)sqrt((x)*(x)+(y)*(y)+(z)*(z));
			}
		}
	}

}

void cal_pc20(float th, float p[][3], float *v, float pc[][3], int *pct)
{
	float pt[2][3];
	float vt[2],pcs[3];
	int pcts;
	int i,j,k;

	for(k=0; k<2; k++){
		for(j=0; j<3; j++){
			for(i=0; i<3; i++){
				pt[0][i] = p[j+(k*4)][i];
		        pt[1][i] = p[(j+1)+(k*4)][i];
			}
	        vt[0] = v[j+(k*4)];
	        vt[1] = v[(j+1)+(k*4)];
	        cal_pc_sub(th, pt, vt, pcs, &pcts);
	        pct[j+(k*4)] = pcts;
	        for(i=0; i<3; i++){
				pc[j+(k*4)][i] = pcs[i];
			}
		}
	}

	for(j=0; j<2; j++){
		for(i=0; i<3; i++){
			pt[0][i] = p[(j*4)][i];
		    pt[1][i] = p[(j*4)+3][i];
		}
	    vt[0] = v[(j*4)];
	    vt[1] = v[(j*4)+3];
	    cal_pc_sub(th, pt, vt, pcs, &pcts);
	    pct[(j*4)+3] = pcts;
	    for(i=0; i<3; i++){
			pc[(j*4)+3][i] = pcs[i];
		}
	}

	for(j=0; j<4; j++){
		for(i=0; i<3; i++){
			pt[0][i] = p[j][i];
		    pt[1][i] = p[j+4][i];
		}
	    vt[0] = v[j];
	    vt[1] = v[j+4];
	    cal_pc_sub(th, pt, vt, pcs, &pcts);
	    pct[j+8] = pcts;
	    for(i=0; i<3; i++){
			pc[j+8][i] = pcs[i];
		}
	}

	for(j=0; j<2; j++){
		for(i=0; i<3; i++){
			pt[0][i] = p[(j*3)+1][i];
		    pt[1][i] = p[(j*3)+3][i];
		}
	    vt[0] = v[(j*3)+1];
	    vt[1] = v[(j*3)+3];
	    cal_pc_sub(th, pt, vt, pcs, &pcts);
	    pct[j+12] = pcts;
	    for(i=0; i<3; i++){
			pc[j+12][i] = pcs[i];
		}
	}


	for(j=0; j<2; j++){
		for(i=0; i<3; i++){
			pt[0][i] = p[1][i];
		    pt[1][i] = p[j*2+4][i];
		}
	    vt[0] = v[1];
	    vt[1] = v[j*2+4];
	    cal_pc_sub(th, pt, vt, pcs, &pcts);
	    pct[j+14] = pcts;
	    for(i=0; i<3; i++){
			pc[j+14][i] = pcs[i];
		}
	}

	for(j=0; j<2; j++){
		for(i=0; i<3; i++){
			pt[0][i] = p[3][i];
		    pt[1][i] = p[6-(j*2)][i];
		}
	    vt[0] = v[3];
	    vt[1] = v[6-(j*2)];
	    cal_pc_sub(th, pt, vt, pcs, &pcts);
	    pct[j+16] = pcts;
	    for(i=0; i<3; i++){
			pc[j+16][i] = pcs[i];
		}
	}

}

void cal_pc21(float th, float p[][3], float *v, float pc[][3], int *pct)
{
	float pt[2][3];
	float vt[2],pcs[3];
	int pcts;
	int i,j,k;

	for(k=0; k<2; k++){
		for(j=0; j<3; j++){
			for(i=0; i<3; i++){
				pt[0][i] = p[j+(k*4)][i];
		        pt[1][i] = p[(j+1)+(k*4)][i];
			}
	        vt[0] = v[j+(k*4)];
	        vt[1] = v[(j+1)+(k*4)];
	        cal_pc_sub(th, pt, vt, pcs, &pcts);
	        pct[j+(k*4)] = pcts;
	        for(i=0; i<3; i++){
				pc[j+(k*4)][i] = pcs[i];
			}
		}
	}

	for(j=0; j<2; j++){
		for(i=0; i<3; i++){
			pt[0][i] = p[(j*4)][i];
		    pt[1][i] = p[(j*4)+3][i];
		}
	    vt[0] = v[(j*4)];
	    vt[1] = v[(j*4)+3];
	    cal_pc_sub(th, pt, vt, pcs, &pcts);
	    pct[(j*4)+3] = pcts;
	    for(i=0; i<3; i++){
			pc[(j*4)+3][i] = pcs[i];
		}
	}

	for(j=0; j<4; j++){
		for(i=0; i<3; i++){
			pt[0][i] = p[j][i];
		    pt[1][i] = p[j+4][i];
		}
	    vt[0] = v[j];
	    vt[1] = v[j+4];
	    cal_pc_sub(th, pt, vt, pcs, &pcts);
	    pct[j+8] = pcts;
	    for(i=0; i<3; i++){
			pc[j+8][i] = pcs[i];
		}
	}

	for(j=0; j<2; j++){
		for(i=0; i<3; i++){
			pt[0][i] = p[(j*5)+0][i];
		    pt[1][i] = p[(j*5)+0][i];
		}
	    vt[0] = v[(j*5)+0];
	    vt[1] = v[(j*5)+2];
	    cal_pc_sub(th, pt, vt, pcs, &pcts);
	    pct[j+12] = pcts;
	    for(i=0; i<3; i++){
			pc[j+12][i] = pcs[i];
		}
	}


	for(j=0; j<2; j++){
		for(i=0; i<3; i++){
			pt[0][i] = p[5][i];
		    pt[1][i] = p[j*2+0][i];
		}
	    vt[0] = v[5];
	    vt[1] = v[j*2+0];
	    cal_pc_sub(th, pt, vt, pcs, &pcts);
	    pct[j+14] = pcts;
	    for(i=0; i<3; i++){
			pc[j+14][i] = pcs[i];
		}
	}

	for(j=0; j<2; j++){
		for(i=0; i<3; i++){
			pt[0][i] = p[7][i];
		    pt[1][i] = p[2-(j*2)][i];
		}
	    vt[0] = v[7];
	    vt[1] = v[2-(j*2)];
	    cal_pc_sub(th, pt, vt, pcs, &pcts);
	    pct[j+16] = pcts;
	    for(i=0; i<3; i++){
			pc[j+16][i] = pcs[i];
		}
	}

}

void count_pc2(int *pct, int *num)
{
	int i;

	(*num) = 0;
	for(i=0; i<6; i++){
		(*num)+=(int)pct[i];
	}
}

void mk_faces2(FILE *fp, float pcs[][3], int *pcts, float *fc, float tr)
{
	float pf[4][3];
	int n;
	int i,j,k;
	
	count_pc2(pcts, &n);


    if(n<3){
		return;
	}else if(n==3){
		k = 0;
		for(i=0; i<6; i++){
			if(pcts[i]==1){
				for(j=0; j<3; j++){
					pf[k][j] = pcs[i][j];
				}
				k++;
			    if(k==3){
					dr_face(fp, pf, n, fc, tr);
					return;
				}
			}
		}
	}else if(n==4){
		if(pcts[0]<1){
			for(i=0; i<3; i++){
				pf[0][i] = pcs[1][i];
			    pf[1][i] = pcs[2][i];
			    pf[2][i] = pcs[5][i];
			    pf[3][i] = pcs[3][i];
			}
			dr_face(fp, pf, n, fc, tr);
			return;
		}else{
			if(pcts[1]<1){
				for(i=0; i<3; i++){
					pf[0][i] = pcs[0][i];
			        pf[1][i] = pcs[3][i];
			        pf[2][i] = pcs[4][i];
			        pf[3][i] = pcs[2][i];
				}
				dr_face(fp, pf, n, fc, tr);
			    return;
			}else{
				for(i=0; i<3; i++){
					pf[0][i] = pcs[0][i];
			        pf[1][i] = pcs[1][i];
			        pf[2][i] = pcs[4][i];
			        pf[3][i] = pcs[5][i];
				}
				dr_face(fp, pf, n, fc, tr);
			    return;
			}
		}
	}else{
		printf("error!");
		getchar();
	}

}

void cal_pc3(float th, float p[][3], float *v, float pc[][3], int *pct)
{
	float pt[2][3];
	float vt[2],pcs[3];
	int pcts;
	int i,j;

	for(j=0; j<3; j++){
		for(i=0; i<3; i++){
			pt[0][i] = p[0][i];
		    pt[1][i] = p[j+1][i];
		}
	    vt[0] = v[0];
	    vt[1] = v[j+1];
	    cal_pc_sub(th, pt, vt, pcs, &pcts);
	    pct[j] = pcts;
	    for(i=0; i<3; i++){
			pc[j][i] = pcs[i];
		}
	}

	for(j=0; j<2; j++){
		for(i=0; i<3; i++){
			pt[0][i] = p[j+1][i];
		    pt[1][i] = p[j+2][i];
		}
	    vt[0] = v[j+1];
	    vt[1] = v[j+2];
	    cal_pc_sub(th, pt, vt, pcs, &pcts);
	    pct[j+3] = pcts;
	    for(i=0; i<3; i++){
			pc[j+3][i] = pcs[i];
		}
	}

	for(i=0; i<3; i++){
			pt[0][i] = p[3][i];
		    pt[1][i] = p[1][i];
	}
	vt[0] = v[3];
	vt[1] = v[1];
	cal_pc_sub(th, pt, vt, pcs, &pcts);
	pct[5] = pcts;
	for(i=0; i<3; i++){
		pc[5][i] = pcs[i];
	}

}

void mk_pt_index(int pind[][6])
{
	pind[0][0] = 0;
	pind[0][1] = 8;
	pind[0][2] = 3;
	pind[0][3] = 14;
	pind[0][4] = 17;
	pind[0][5] = 12;

	pind[1][0] = 5;
	pind[1][1] = 4;
	pind[1][2] = 9;
	pind[1][3] = 13;
	pind[1][4] = 14;
	pind[1][5] = 15;

	pind[2][0] = 1;
	pind[2][1] = 2;
	pind[2][2] = 10;
	pind[2][3] = 12;
	pind[2][4] = 16;
	pind[2][5] = 15;

	pind[3][0] = 6;
	pind[3][1] = 11;
	pind[3][2] = 7;
	pind[3][3] = 16;
	pind[3][4] = 17;
	pind[3][5] = 13;

	pind[4][0] = 12;
	pind[4][1] = 15;
	pind[4][2] = 14;
	pind[4][3] = 16;
	pind[4][4] = 13;
	pind[4][5] = 17;
}

void mk_equ_face2_sub(FILE *fp, float th, float p[][3], float *v, float *fc, float tr)
{
	float pc[18][3];
	int pct[18];
	float pctt[6][3];
	int pcttt[6];
	int pind[5][6];
	int i,j,k;

	cal_pc20(th, p, v, pc, pct);

	mk_pt_index(pind);

	for(i=0; i<5; i++){
		for(j=0; j<6; j++){
			for(k=0; k<3; k++){
				pctt[j][k] = pc[(pind[i][j])][k];
			}
			pcttt[j] = pct[(pind[i][j])];
		}
		mk_faces2(fp, pctt, pcttt, fc, tr);
	}

}

void mk_equ_face2(FILE *fp, float th, float v[][nY][nZ], float *fc, float tr)
{
	float x,y,z;
	float vt[8];
	float p[8][3];

	for(x=0; x<(NX)*D; x+=D){
		for(y=0; y<(NY)*D; y+=D){
			for(z=0; z<(NZ)*D; z+=D){
				mk_cube(x, y, z, p, v, vt);
				mk_equ_face2_sub(fp, th, p, vt, fc, tr);
			}
		}
	}
}

void clock_start(clock_t *start)
{
	*start = clock();
}

void clock_end(clock_t start)
{
	double time;

	time = (clock()-start)/(double)CLOCKS_PER_SEC;
	printf("time = %lfs\n", time);

}

void trans_v(float *v, float *vtr, float p[][3], float ptr[][3])
{
	int i;

	vtr[0] = v[1];
	vtr[1] = v[2];
	vtr[2] = v[3];
	vtr[3] = v[0];
	vtr[4] = v[5];
	vtr[5] = v[6];
	vtr[6] = v[7];
	vtr[7] = v[4];

	for(i=0; i<3; i++){
	    ptr[0][i] = p[1][i];
	    ptr[1][i] = p[2][i];
	    ptr[2][i] = p[3][i];
	    ptr[3][i] = p[0][i];
	    ptr[4][i] = p[5][i];
	    ptr[5][i] = p[6][i];
	    ptr[6][i] = p[7][i];
	    ptr[7][i] = p[4][i];
	}
}


void mk_cube3(float x, float y, float z, float dx, float dy, float dz, float p[][3], float v[][3][3], float *vt)
{
	int i;
	int xt,yt,zt;

	for(i=0; i<2; i++){
		p[0+(i*4)][0] = dx;
		p[0+(i*4)][1] = dy;
		p[0+(i*4)][2] = dz+i*D;

		p[1+(i*4)][0] = dx+D;
		p[1+(i*4)][1] = dy;
		p[1+(i*4)][2] = dz+i*D;

		p[2+(i*4)][0] = dx+D;
		p[2+(i*4)][1] = dy+D;
		p[2+(i*4)][2] = dz+i*D;

		p[3+(i*4)][0] = dx;
		p[3+(i*4)][1] = dy+D;
		p[3+(i*4)][2] = dz+i*D;
	}
	for(i=0; i<8; i++){
		xt = (int)(p[i][0]-x)/(int)D;
		yt = (int)(p[i][1]-y)/(int)D;
		zt = (int)(p[i][2]-z)/(int)D;
		vt[i] = v[xt][yt][zt];
	}

}

void mk_cube2(float x, float y, float z, float v[][nY][nZ], float vt[][3][3])
{
	int i,j,k;
	int xt,yt,zt;

	for(i=0; i<3; i++){
		for(j=0; j<3; j++){
			for(k=0; k<3; k++){
				xt = i+(int)(x/D);
				yt = j+(int)(y/D);
				zt = k+(int)(z/D);
				vt[i][j][k] = v[xt][yt][zt];
			}
		}
	}

}


void mk_equ_face001(FILE *fp, float th, float v[][nY][nZ], float *fc, float tr)
{
	float x,y,z;
	float vt[8];
	float p[8][3];
	int rm;

	for(x=0; x<(NX)*D; x+=D){
		for(y=0; y<(NY)*D; y+=D){
			for(z=0; z<(NZ)*D; z+=D){
				mk_cube(x, y, z, p, v, vt);
				rm_cube(th, vt, &rm);
				if(rm == 1 ) continue;
				mk_equ_face_sub(fp, th, p, vt, fc, tr);
			}
		}
	}
}

void rm_cube(float th, float *vt, int *rm)
{
	int i,uc,dc;

	uc = 0;
	dc = 0;
	for(i=0; i<8; i++){
		if( vt[i]>=th ){
			uc++;
		}else if( vt[i]<th ){
			dc++;
		}
	}

	if( (uc==8) || (dc==8) ){
		(*rm) = 1;
	}else{
		(*rm) = 0;
	}
}

void mk_equ_face0021(FILE *fp, float x, float y, float z, float th, float v[][nY][nZ], float *fc, float tr)
{
	float dx,dy,dz;
	int xt,yt,zt,nt;
	float vt[3][3][3],vtt[8],vtr[8];
	float p[8][3],ptr[8][3];
	int rm;


	mk_cube2(x, y, z, v, vt);

	for(dx=x; dx<(x+2*D); dx+=D){
		for(dy=y; dy<(y+2*D); dy+=D){
			for(dz=z; dz<(z+2*D); dz+=D){
				mk_cube3(x, y, z, dx, dy, dz, p, vt, vtt);


				rm_cube(th, vtt, &rm);
				if( rm == 1 ) continue;

		        xt = (int)((dx-x)/D);
		        yt = (int)((dy-y)/D);
		        zt = (int)((dz-z)/D);
				nt = xt+yt+zt;
				if( (nt%2)==1 ){
					trans_v(vtt, vtr, p, ptr);
					mk_equ_face2_sub(fp, th, ptr, vtr, fc, tr);
				}else{
					mk_equ_face2_sub(fp, th, p, vtt, fc, tr);
				}
			}
		}
	}
}

void mk_equ_face002(FILE *fp, float th, float v[][nY][nZ], float *fc, float tr)
{
	float x,y,z;

	for(x=0; x<(XM-D); x+=(2*D)){
		for(y=0; y<(YM-D); y+=(2*D)){
			for(z=0; z<(ZM-D); z+=(2*D)){
				mk_equ_face0021(fp, x, y, z, th, v, fc, tr);
			}
		}
	}
}